1 Riassunti Numerici per Variabili Quantitative

Quando si dispone dei dati per una variabile quantitativa è necessario sintetizzarli in modo da far emergere i fatti salienti, rinunciando ai dettagli relativi alle singole unità per dare una lettura di insieme del fenomeno.

I due aspetti su cui si concentra l’attenzione sono:

  • la tendenza centrale: con un singolo valore posso riassumere l’ordine di grandezza del fenomeno?
  • la dispersione (o variabilità): con un singolo valore posso fornire un’idea di quanto i dati siano diversi tra loro? (o simili)

Nel seguito useremo in particolare la variabile HEIGHT del dataset Esempio.txt.

dati <- read.table("Esempio.txt", header = TRUE, sep = ",")
head(dati)
str(dati)
#> 'data.frame':    92 obs. of  4 variables:
#>  $ HEIGHT  : num  66 72 73.5 73 69 73 72 74 72 71 ...
#>  $ WEIGHT  : int  140 145 160 190 155 165 150 190 195 138 ...
#>  $ ACTIVITY: int  2 2 3 1 2 1 3 2 2 2 ...
#>  $ SMOKES  : int  2 2 1 1 2 2 2 2 2 2 ...
x <- dati$HEIGHT
n    <- length(x)

⚠️ La variabile HEIGHT non presenta NA. In presenza di valori mancanti è sempre necessario aggiungere na.rm = TRUE nelle funzioni statistiche.

2 Indici di Tendenza Centrale

Gli indici di tendenza centrale sono anche noti come medie, distinte in medie analitiche e medie di posizione.

2.1 Medie Analitiche

2.1.1 Media Aritmetica

\[M = \sum_{i=1}^{n} \frac{x_i}{n}\]

Proprietà principali:

  • assume un valore compreso fra \(\min(x_i)\) e \(\max(x_i)\)
  • la somma degli scarti dalla media è nulla: \(\sum_{i=1}^{n}(x_i - M) = 0\), da cui \(\sum x_i = nM\) (la media è il valore che sostituito a ogni dato lascia invariata la somma)
  • minimizza la somma degli scarti quadratici: \(M = \arg\min_{a \in \mathbb{R}} \sum_{i=1}^{n}(x_i - a)^2\)
M <- mean(x)
M
#> [1] 68.71739
# Verifica proprietà a: min <= M <= max
M >= min(x) & M <= max(x)
#> [1] TRUE
# Verifica proprietà b: scarti nulli
round(sum(x - M), 10)
#> [1] 0
# Verifica proprietà c: M minimizza la somma degli scarti quadratici
a       <- sort(c(seq(min(x), max(x), 0.5), mean(x)))
varTest <- sapply(a, function(a) mean((x - a)^2))
min.a   <- which.min(varTest)

plot(a, varTest, type = "b",
     xlab = "a", ylab = expression(sum((x[i]-a)^2)),
     main = "La media minimizza la somma degli scarti quadratici")
points(a[min.a], varTest[min.a], pch = 20, col = "red", cex = 1.5)

💡 In R, il comando plot() permette di ottenere un diagramma Cartesiano; in questo caso riportiamo sull’asse x i valori di a e sull’asse y i valori di varTest

💡 In R, il comando points() permette usando il sistema di coordinate (x,y) di aggiungere un punto al grafico

2.1.2 Altre Medie Analitiche

Partendo dalla stessa logica — cercare il valore che lascia invariata una certa funzione dei dati — si ottengono:

Media Formula Invariante
Geometrica \(M_g = \left[\prod x_i\right]^{1/n}\) prodotto dei dati
Armonica \(M_{ar} = n / \sum(1/x_i)\) somma dei reciproci
Di potenza \(r\) \(M_r = \left[\sum x_i^r / n\right]^{1/r}\) somma delle potenze \(r\)-esime

Proprietà:

  • \(M_{ar} \leq M_g \leq M \leq M_2\) (uguaglianza solo se tutti i dati sono uguali).
  • Vale che \(M^2_2 \geq M^2\)
  • La potenza \(r\)-esima di \(M_r\) fornisce il momento empirico \(r\)-esimo.
# Media geometrica
Mg <- prod(x)^(1/n)
Mg
#> [1] 68.62003
# Media armonica
Ma <- n / sum(1/x)
Ma
#> [1] 68.52173
# Media di potenza (r = 2)
r  <- 2
M2 <- (sum(x^r) / n)^(1/r)
M2
#> [1] 68.8137
# Verifica ordinamento
Ma <= Mg & Mg <= M & M <= M2
#> [1] TRUE
M2^2 >= M^2
#> [1] TRUE

2.2 Quantili e Medie di Posizione

Fissata una proporzione \(p \in [0,1]\), il quantile empirico di ordine \(p\) è quel valore \(x_p \in \mathbb{R}\) tale che:

\[x_p : \frac{\#(x_i \leq x_p)}{n} = p\]

Considerando i dati ordinati \(x_{(1)} \leq \ldots \leq x_{(n)}\), il valore \(x_{(i)}\) rappresenta il quantile di ordine \(p = i/n\).

💡 In R quantile() usa l’interpolazione lineare fra i due quantili osservati più vicini. La convenzione è che il minimo sia il quantile di ordine 0 e il massimo quello di ordine 1. Si veda ?quantile per le diverse strategie disponibili.

Altre convenzioni: * Se disponiamo di un insieme di dati non particolarmente numeroso potremmo chiedere di contare invece dei dati minori o uguali, quelli strettamente minori di \(x_p\): la sequenza ordinata fornirebbe i quantili \(\frac{0}{n}, \frac{1}{n}, \dots , \frac{n-2}{n} , \frac{n-1}{n}\).

  • A volte si usa una diversa convenzione per cui i dati ordinati \(x_{(1)} \leq x_{(2)} \leq \dots, \leq x_{(n-1)} \leq x_{(n)}\) rappresentano i quantili di ordine \(\frac{0.5}{n}, \frac{1.5}{n}, \dots , \frac{n-1.5}{n} , \frac{n-0.5}{n}\). Cioè secondo tale convenzione \(x_{(i)}\) rappresenta il quantile di ordine \(x_{p}\), con \(p=\frac{i-0.5}{n}\).

  • Se \(n\) è grande ci si attende che le differenze derivanti dalle diverse convenzioni si annullino.

# Le tre convenzioni principali per i quantili empirici
cbind(sort(x),
      CONV1 = 1:n / n,
      CONV2 = 0:(n-1) / n,
      CONV3 = (0.5:(n-0.5)) / n) 
#>                  CONV1      CONV2       CONV3
#>  [1,] 61.00 0.01086957 0.00000000 0.005434783
#>  [2,] 61.75 0.02173913 0.01086957 0.016304348
#>  [3,] 62.00 0.03260870 0.02173913 0.027173913
#>  [4,] 62.00 0.04347826 0.03260870 0.038043478
#>  [5,] 62.00 0.05434783 0.04347826 0.048913043
#>  [6,] 62.00 0.06521739 0.05434783 0.059782609
#>  [7,] 62.75 0.07608696 0.06521739 0.070652174
#>  [8,] 63.00 0.08695652 0.07608696 0.081521739
#>  [9,] 63.00 0.09782609 0.08695652 0.092391304
#> [10,] 63.00 0.10869565 0.09782609 0.103260870
#> [11,] 63.00 0.11956522 0.10869565 0.114130435
#> [12,] 64.00 0.13043478 0.11956522 0.125000000
#> [13,] 64.00 0.14130435 0.13043478 0.135869565
#> [14,] 65.00 0.15217391 0.14130435 0.146739130
#> [15,] 65.00 0.16304348 0.15217391 0.157608696
#> [16,] 65.00 0.17391304 0.16304348 0.168478261
#> [17,] 65.00 0.18478261 0.17391304 0.179347826
#> [18,] 65.50 0.19565217 0.18478261 0.190217391
#> [19,] 66.00 0.20652174 0.19565217 0.201086957
#> [20,] 66.00 0.21739130 0.20652174 0.211956522
#> [21,] 66.00 0.22826087 0.21739130 0.222826087
#> [22,] 66.00 0.23913043 0.22826087 0.233695652
#> [23,] 66.00 0.25000000 0.23913043 0.244565217
#> [24,] 66.00 0.26086957 0.25000000 0.255434783
#> [25,] 66.00 0.27173913 0.26086957 0.266304348
#> [26,] 66.00 0.28260870 0.27173913 0.277173913
#> [27,] 67.00 0.29347826 0.28260870 0.288043478
#> [28,] 67.00 0.30434783 0.29347826 0.298913043
#> [29,] 67.00 0.31521739 0.30434783 0.309782609
#> [30,] 67.00 0.32608696 0.31521739 0.320652174
#> [31,] 67.00 0.33695652 0.32608696 0.331521739
#> [32,] 67.00 0.34782609 0.33695652 0.342391304
#> [33,] 67.00 0.35869565 0.34782609 0.353260870
#> [34,] 68.00 0.36956522 0.35869565 0.364130435
#> [35,] 68.00 0.38043478 0.36956522 0.375000000
#> [36,] 68.00 0.39130435 0.38043478 0.385869565
#> [37,] 68.00 0.40217391 0.39130435 0.396739130
#> [38,] 68.00 0.41304348 0.40217391 0.407608696
#> [39,] 68.00 0.42391304 0.41304348 0.418478261
#> [40,] 68.00 0.43478261 0.42391304 0.429347826
#> [41,] 68.00 0.44565217 0.43478261 0.440217391
#> [42,] 68.00 0.45652174 0.44565217 0.451086957
#> [43,] 68.00 0.46739130 0.45652174 0.461956522
#> [44,] 69.00 0.47826087 0.46739130 0.472826087
#> [45,] 69.00 0.48913043 0.47826087 0.483695652
#> [46,] 69.00 0.50000000 0.48913043 0.494565217
#> [47,] 69.00 0.51086957 0.50000000 0.505434783
#> [48,] 69.00 0.52173913 0.51086957 0.516304348
#> [49,] 69.00 0.53260870 0.52173913 0.527173913
#> [50,] 69.00 0.54347826 0.53260870 0.538043478
#> [51,] 69.00 0.55434783 0.54347826 0.548913043
#> [52,] 69.00 0.56521739 0.55434783 0.559782609
#> [53,] 69.00 0.57608696 0.56521739 0.570652174
#> [54,] 69.50 0.58695652 0.57608696 0.581521739
#> [55,] 70.00 0.59782609 0.58695652 0.592391304
#> [56,] 70.00 0.60869565 0.59782609 0.603260870
#> [57,] 70.00 0.61956522 0.60869565 0.614130435
#> [58,] 70.00 0.63043478 0.61956522 0.625000000
#> [59,] 70.00 0.64130435 0.63043478 0.635869565
#> [60,] 70.00 0.65217391 0.64130435 0.646739130
#> [61,] 71.00 0.66304348 0.65217391 0.657608696
#> [62,] 71.00 0.67391304 0.66304348 0.668478261
#> [63,] 71.00 0.68478261 0.67391304 0.679347826
#> [64,] 71.00 0.69565217 0.68478261 0.690217391
#> [65,] 71.00 0.70652174 0.69565217 0.701086957
#> [66,] 71.00 0.71739130 0.70652174 0.711956522
#> [67,] 71.50 0.72826087 0.71739130 0.722826087
#> [68,] 72.00 0.73913043 0.72826087 0.733695652
#> [69,] 72.00 0.75000000 0.73913043 0.744565217
#> [70,] 72.00 0.76086957 0.75000000 0.755434783
#> [71,] 72.00 0.77173913 0.76086957 0.766304348
#> [72,] 72.00 0.78260870 0.77173913 0.777173913
#> [73,] 72.00 0.79347826 0.78260870 0.788043478
#> [74,] 72.00 0.80434783 0.79347826 0.798913043
#> [75,] 72.00 0.81521739 0.80434783 0.809782609
#> [76,] 73.00 0.82608696 0.81521739 0.820652174
#> [77,] 73.00 0.83695652 0.82608696 0.831521739
#> [78,] 73.00 0.84782609 0.83695652 0.842391304
#> [79,] 73.00 0.85869565 0.84782609 0.853260870
#> [80,] 73.00 0.86956522 0.85869565 0.864130435
#> [81,] 73.00 0.88043478 0.86956522 0.875000000
#> [82,] 73.00 0.89130435 0.88043478 0.885869565
#> [83,] 73.50 0.90217391 0.89130435 0.896739130
#> [84,] 73.50 0.91304348 0.90217391 0.907608696
#> [85,] 74.00 0.92391304 0.91304348 0.918478261
#> [86,] 74.00 0.93478261 0.92391304 0.929347826
#> [87,] 74.00 0.94565217 0.93478261 0.940217391
#> [88,] 74.00 0.95652174 0.94565217 0.951086957
#> [89,] 74.00 0.96739130 0.95652174 0.961956522
#> [90,] 75.00 0.97826087 0.96739130 0.972826087
#> [91,] 75.00 0.98913043 0.97826087 0.983695652
#> [92,] 75.00 1.00000000 0.98913043 0.994565217

2.2.1 La Mediana

La mediana \(Me\) è il valore per cui il 50% dei dati è inferiore (o uguale) a esso. Per convenzione:

\[Me = \begin{cases} x_{\left(\frac{n+1}{2}\right)} & n \text{ dispari} \\ \dfrac{x_{\left(\frac{n}{2}\right)} + x_{\left(\frac{n}{2}+1\right)}}{2} & n \text{ pari} \end{cases}\]

La mediana minimizza la somma degli scarti assoluti: \(Me = \arg\min_{a \in \mathbb{R}} \sum_{i=1}^n |x_i - a|\). Per questo è molto meno sensibile ai valori estremi rispetto alla media — è più resistente.

median(x)
#> [1] 69
quantile(x, probs = 0.5)
#> 50% 
#>  69
# Calcolo manuale
xs <- sort(x)
ifelse(n %% 2 == 0, (xs[n/2] + xs[n/2+1])/2, xs[(n+1)/2])
#> [1] 69
# Verifica proprietà di minimo
a       <- sort(c(seq(min(x), max(x), 0.5), median(x)))
AbsTest <- sapply(a, function(a) sum(abs(x - a)))
min.a   <- which.min(AbsTest)
plot(a, AbsTest, type = "b",
     xlab = "a", ylab = expression(sum(abs(x[i]-a))),
     main = "La mediana minimizza la somma degli scarti assoluti")
points(a[min.a], AbsTest[min.a], pch = 20, col = "red", cex = 1.5)

2.2.2 Quartili, Decili e Quantili

I quantili \(x_{0.25}\), \(x_{0.5}\), \(x_{0.75}\) sono detti primo, secondo e terzo quartile. Il 50% dei valori centrali è compreso fra Q1 e Q3.

Q1 <- quantile(x, probs = 0.25); Q1
#> 25% 
#>  66
Q3 <- quantile(x, probs = 0.75); Q3
#> 75% 
#>  72
# Decili
quantile(x, probs = seq(0, 1, 0.1))
#>    0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100% 
#> 61.00 63.00 66.00 67.00 68.00 69.00 70.00 71.00 72.00 73.45 75.00
# Quantili specifici (es. 2.5° e 97.5° percentile)
quantile(x, probs = c(0.025, 0.975))
#>   2.5%  97.5% 
#> 62.000 74.725

2.2.3 Riassunto dei Cinque Valori

\(\langle \min,\, Q_1,\, Me,\, Q_3,\, \max \rangle\) — il five numbers summary — fornisce una panoramica completa sulla distribuzione. summary() aggiunge anche la media.

fivenum(x)
#> [1] 61 66 69 72 75
summary(x)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   61.00   66.00   69.00   68.72   72.00   75.00

⚠️ fivenum() e summary() rimuovono automaticamente gli NA, a differenza delle funzioni statistiche come mean() e var().

2.2.4 Media Sfrondata

La media sfrondata (trimmed mean) esclude una proporzione \(\alpha\) di valori estremi per lato prima di calcolare la media, rendendola meno sensibile agli outlier.

mean(x, trim = 0.1)    # esclude il 10% più basso e il 10% più alto
#> [1] 68.83784
mean(sort(x)[10:83])
#> [1] 68.83784
mean(x)                # confronto con la media ordinaria
#> [1] 68.71739

3 Indici di Dispersione

3.1 Varianza e Scarto Quadratico Medio

La devianza è la somma degli scarti quadratici dalla media: \(DEV = \sum_{i=1}^n (x_i - M)^2\).

La varianza è la devianza media: \(V = DEV/n\).

⚠️ In R var() divide per \(n-1\) (capiremo il perchè duranto il corso di inferenza statistica). Per \(n\) grande la differenza è trascurabile, ma è bene esserne consapevoli.

Lo scarto quadratico medio (SQM) o deviazione standard \(sd = \sqrt{V}\) riporta l’indice nella stessa unità di misura della variabile.

var(x)
#> [1] 13.39041
sqrt(var(x))
#> [1] 3.659291
sd(x)
#> [1] 3.659291
# Varianza e sd per fumatori (SMOKES==1) e non fumatori (SMOKES==2)
var(x[dati$SMOKES == 1]);  sd(x[dati$SMOKES == 1])
#> [1] 16.57564
#> [1] 4.07132
var(x[dati$SMOKES == 2]);  sd(x[dati$SMOKES == 2])
#> [1] 12.1767
#> [1] 3.489512
# Media per confronto
mean(x[dati$SMOKES == 1])
#> [1] 69.02679
mean(x[dati$SMOKES == 2])
#> [1] 68.58203

📌 I due gruppi hanno media simile ma i fumatori mostrano varianza maggiore. Poiché media e varianza non sono resistenti, è utile affiancare misure robuste.

3.2 Scarto Interquartile, MAD e Coefficiente di Variazione

Indice Formula Resistente? Funzione R
Scarto interquartile \(IQR = Q_3 - Q_1\) Sì IQR()
MAD \(\text{Mediana}(|x_i - Me|)\) Sì mad(constant=1)
Coeff. di variazione \(CV = sd / M\) No — 
SI  <- c(IQR(x[dati$SMOKES == 1]),
         IQR(x[dati$SMOKES == 2]))

MAD <- c(mad(x[dati$SMOKES == 1], constant = 1),
         mad(x[dati$SMOKES == 2], constant = 1))
mad(x, constant = 1)
#> [1] 3
1 * median(abs(x-median(x)))
#> [1] 3
CV  <- c(sd(x[dati$SMOKES == 1]) / mean(x[dati$SMOKES == 1]),
         sd(x[dati$SMOKES == 2]) / mean(x[dati$SMOKES == 2]))

rbind(SI, MAD, CV)
#>           [,1]       [,2]
#> SI  6.25000000 5.25000000
#> MAD 4.00000000 3.00000000
#> CV  0.05898174 0.05088085

📌 Lo scarto interquartile e il MAD sono simili nei due gruppi, indicando che la differenza di varianza tra fumatori e non fumatori non è dovuta esclusivamente a valori anomali.

4 Simmetria e Curtosi

4.1 Asimmetria

Una distribuzione presenta asimmetria positiva (coda destra lunga) se i dati a destra della mediana sono più dispersi; asimmetria negativa (coda sinistra lunga) nel caso opposto. I decili permettono di valutare visivamente l’asimmetria confrontando le distanze dalla mediana verso destra e verso sinistra.

quantile(x, probs = seq(0, 1, 0.1), digits = 2)
#>    0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100% 
#> 61.00 63.00 66.00 67.00 68.00 69.00 70.00 71.00 72.00 73.45 75.00
quantile(dati$WEIGHT, probs = seq(0, 1, 0.1), digits = 2)
#>    0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100% 
#>  95.0 116.0 123.4 130.0 138.8 145.0 150.0 155.0 160.0 179.5 215.0

📌 Per HEIGHT: distanza mediana–minimo = 8, mediana–massimo = 6 → leggera asimmetria sinistra. Distanza mediana - primo decile è 6, mediana - nono decile è 4.45 \(\implies\) Coda sinistra leggermente più lunga - Distribuzione leggermente asimmetrica a sinistra (ma molto vicina alla simmetria)

Per WEIGHT: distanza mediana–minimo = 50, mediana–massimo = 70; → leggera asimmetria destra; distanza tra mediana - primo decile è 29, quella tra mediana - nono decile è 34.5 \(\implies\) Coda destra leggermente più lunga - Distribuzione asimmetrica a destra

4.1.1 Indice di Galton e Indice \(K\)

Due indici di simmetria sono \[G = \frac{Q_3 + Q_1 - 2Q_2}{Q_3 - Q_1}\]

\[K = \frac{x_{0.9} + x_{0.1} - 2x_{0.5}}{x_{0.9} - x_{0.1}}\]

\(G\) cattura l’asimmetria nella parte centrale; \(K\) include informazione sulle code. Entrambi sono nulli in caso di simmetria, positivi per asimmetria destra e negativi per asimmetria sinistra.

# Indice di Galton
Qh <- fivenum(x, na.rm = TRUE)
Qw <- fivenum(dati$WEIGHT, na.rm = TRUE)

G_HEIGHT <- (Qh[4] + Qh[2] - 2*Qh[3]) / (Qh[4] - Qh[2]); G_HEIGHT
#> [1] 0
G_WEIGHT <- (Qw[4] + Qw[2] - 2*Qw[3]) / (Qw[4] - Qw[2]); G_WEIGHT
#> [1] -0.2903226
# Indice K (basato su decili)
Dh <- quantile(x,           probs = seq(0, 1, 0.1))
Dw <- quantile(dati$WEIGHT, probs = seq(0, 1, 0.1))

K_HEIGHT <- (Dh[10] + Dh[2] - 2*Dh[6]) / (Dh[10] - Dh[2])
K_HEIGHT
#>        90% 
#> -0.1483254
K_WEIGHT <- (Dw[10] + Dw[2] - 2*Dw[6]) / (Dw[10] - Dw[2])
K_WEIGHT
#>        90% 
#> 0.08661417

Per altezza \(G\) fornisce una indicazione di simmetria, mentre \(K\) fornisce una indicazione di leggera asimmetria sinistra

Per peso, \(G\) fornisce una indicazione di asimmetria a sinistra (se guardo ai dati centrali della distribuzione) mentre fornisce una indicazione (debole) di asimmetria destra (se includo l’informazione sulle code)

4.1.2 Indice di Pearson e Indice di Asimmetria \(\gamma\)

Il coefficiente di Pearson \(\frac{3(M - Me)}{sd}\) sfrutta la distanza tra media e mediana: se la media è superiore alla mediana c’è asimmetria positiva.

L’indice di asimmetria \(\gamma\) è basato sugli scarti dalla media al cubo:

\[\gamma = \frac{\sum_{i=1}^n (x_i - M)^3 / n}{sd^3}\]

# Pearson
3*(mean(x) - median(x)) / sd(x)
#> [1] -0.2316914
3*(mean(dati$WEIGHT) - median(dati$WEIGHT)) / sd(dati$WEIGHT)
#> [1] 0.01923055
# Indice gamma (pacchetto moments)
library(moments)
skewness(x)
#> [1] -0.2209687
skewness(dati$WEIGHT)
#> [1] 0.3639218

📌 Per HEIGHT: leggera asimmetria negativa (confermata da tutti gli indici). Per WEIGHT: asimmetria positiva moderata — la coda destra è più lunga.

4.2 Curtosi

La curtosi misura la pesantezza delle code rispetto alla distribuzione Normale (riferimento: \(\delta = 3\)):

\[\delta = \frac{\sum_{i=1}^n (x_i - M)^4 / n}{sd^4}\]

  • \(\delta < 3\): code più corte della gaussiana (leptocurtosi)
  • \(\delta > 3\): code più pesanti (platicurtosi)
kurtosis(x, na.rm = TRUE)
#> [1] 2.17413

📌 Sebbene sia stata evidenziata una leggera asimmetria a sinistra, l’analisi della curtosi evidenzia che la distribuzione ha code più leggere rispetto a quelle della normale.

5 Moda

La moda è l’unico indice di tendenza centrale applicabile a variabili categoriali. Si definisce diversamente a seconda del tipo di variabile:

  • variabile categoriale: modalità con frequenza relativa massima
  • variabile discreta: valore con frequenza massima
  • variabile continua: si raggruppano i dati in classi e si individua la classe modale (quella con densità di frequenza massima, non frequenza relativa massima)
# Moda per variabile categoriale
table(dati$ACTIVITY) # Moda: modalità 2
#> 
#>  0  1  2  3 
#>  1  9 61 21
# Moda per variabile quantitativa (intervalli di uguale ampiezza)
x_class <- cut(x, breaks = 5, include.lowest = T)
table(x_class) # Moda: (66.6,69.4] 
#> x_class
#>   [61,63.8] (63.8,66.6] (66.6,69.4] (69.4,72.2]   (72.2,75] 
#>          11          15          27          22          17

💡 In distribuzioni unimodali la posizione relativa di moda, mediana e media indica l’asimmetria: con asimmetria positiva (coda destra) \(Mo < Me < M\); con asimmetria negativa l’ordine si inverte.

Se una distribuzione mostra più di un punto di massimo si parla di distribuzione multimodale (bimodale se i massimi sono due). Come esempio classico:

library(MASS)
data(geyser)
hist(geyser$waiting,
     breaks  = c(43, seq(48, 108, by = 5)),
     prob    = TRUE,
     main    = "Tempi di attesa fra eruzioni (geyser)",
     xlab    = "Minuti", col = "lightblue", border = "white")

6 Stima della Densità con il Metodo del Nucleo

L’istogramma dipende dalla scelta delle classi. Il metodo del nucleo (kernel density estimation) calcola la densità in ogni punto \(x\) direttamente dai dati:

\[d(x) = \frac{\#\text{ valori in } (x - h/2,\; x + h/2]}{hn}\] Idea: la densità in \(x\) è proporzionale al numero di osservazioni vicine a \(x\), dove “vicine” è definito da un intorno di ampiezza \(h\)

6.1 Nucleo Rettangolare

Equivale a posizionare un rettangolo di ampiezza \(h\) centrato su ogni \(x_i\) e sommare i contributi:

\[W(x)=\begin{cases} 1 & \mbox{se}~~~ |x| < \frac{1}{2} \\ 0 & \mbox{altrimenti} \end{cases}\]

\[d(x) = \frac{1}{nh} \sum_{i=1}^{n} W\!\left(\frac{x - x_i}{h}\right)\]

# Utilizziamo questi nove dati di esempio 
xx <- c(148, 155, 162, 168, 172, 174, 175, 178, 186)

rett <- function(x, xx, h = 1) {
  n <- length(xx)
  r <- 0
  for (i in 1:n) r <- r + (abs(x - xx[i]) <= h/2) / (n*h)
  r
}

xr <- seq(135, 200, 0.5)

plot(xx, rep(0, 9), ylim = c(0, 0.065), xlim = c(135, 200),
     pch = 4, cex = 1.5, main = "Nucleo rettangolare (h = 10)",
     xlab = "x", ylab = "densità")
abline(h = 0)
for (i in seq_along(xx)) {
  yc <- jitter(1/180)
  segments(xx[i]-10, yc, xx[i]+10, yc, col = i+1, lwd = 3)
  segments(xx[i]-10, 0,  xx[i]-10, yc, col = i+1, lwd = 2, lty = 2)
  segments(xx[i]+10, 0,  xx[i]+10, yc, col = i+1, lwd = 2, lty = 2)
}
points(xr, rett(xr, xx, 20), lwd = 2.5, type = "s")

6.2 Nucleo Gaussiano

Il nucleo gaussiano produce una curva continua e derivabile: ogni \(x_i\) contribuisce con una piccola gaussiana centrata su di esso.

\[d(x) = \frac{1}{nh} \sum_{i=1}^{n} K\!\left(\frac{x - x_i}{h}\right), \quad K(u) = \frac{1}{\sqrt{2\pi}} e^{-u^2/2}\]

nuclg <- function(x, xx, h = 1)
  sapply(x, function(xi) sum(dnorm(xi, xx, h) / length(xx)))

plot(xx, rep(0, 9), ylim = c(0, 0.065), xlim = c(135, 200),
     pch = 4, cex = 1.5, main = "Nucleo gaussiano (h = 4)",
     xlab = "x", ylab = "densità")
abline(h = 0)
for (i in seq_along(xx))
  curve(dnorm(x, xx[i], 4) / length(xx), col = i+1, lwd = 2, lty = 2, add = TRUE)
lines(xr, nuclg(xr, xx, 4), lwd = 2.5)

6.3 Il Parametro di Lisciamento bw

  • \(h\) piccolo → curva irregolare, cattura variazioni locali
  • \(h\) grande → curva liscia ma può nascondere caratteristiche reali (es. bimodalità)
par(mfrow = c(1, 3))
for (h in c(2.5, 4, 6)) {
  plot(xx, rep(0, 9), ylim = c(0, 0.065), xlim = c(135, 200),
       pch = 4, cex = 1.5, main = paste("h =", h), xlab = "x", ylab = "densità")
  abline(h = 0)
  for (i in seq_along(xx))
    curve(dnorm(x, xx[i], h)/length(xx), col = i+1, lwd = 2, lty = 2, add = TRUE)
  lines(xr, nuclg(xr, xx, h), lwd = 2.5)
}

par(mfrow = c(1, 1))

6.4 La Funzione density() in R

density() implementa il nucleo gaussiano con bandwidth ottimale automatica: \[h = 0.9 \cdot A \cdot n^{-1/5}, \quad A = \min\!\left(sd,\; \frac{IQR}{1.34}\right)\]

💡 La scelta del nucleo è meno cruciale di quella di bw: nuclei diversi producono curve simili, mentre bw influenza significativamente il risultato.

den <- density(geyser$waiting)
plot(den, main = paste("Tempi di attesa eruzioni — bw =", round(den$bw, 2)))

par(mfrow = c(1, 3))
plot(density(geyser$waiting, bw =  1), main = "bw = 1 (troppo piccolo)")
plot(density(geyser$waiting, bw =  4), main = "bw = 4 (adeguato)")
plot(density(geyser$waiting, bw = 10), main = "bw = 10 (troppo grande)")

par(mfrow = c(1, 1))

⚠️ Con bw troppo grande si perde la bimodalità dei tempi di attesa. Con bw troppo piccolo la curva è eccessivamente irregolare. Confrontare visivamente più valori prima di scegliere.

7 Trasformazione delle Variabili

7.1 Trasformazioni Lineari e Standardizzazione

Data la trasformazione \(Z_i = bX_i + a\):

  • Media della variabile \(Z\): \(\bar Z = b \bar X + a\)
  • Varianza della variabile \(Z\): \(\sigma^2_Z = b^2 \sigma^2_X\)

La standardizzazione \(Z_i = (X_i - \bar X) / \sigma_X\) produce una variabile con media 0 e varianza 1, utile per confrontare variabili su scale diverse.

z <- (x - mean(x)) / sd(x)
round(mean(z), 10)   # → 0
#> [1] 0
round(var(z),  10)   # → 1
#> [1] 1

7.2 Trasformazioni Non Lineari

Quando una variabile ha code molto pesanti o distribuzione fortemente asimmetrica, la trasformazione logaritmica \(Z_i = \log(X_i)\) rende la distribuzione più simmetrica e la rappresentazione grafica più leggibile.

⚠️ La media della variabile trasformata non è uguale alla trasformazione della media originale (disuguaglianza di Jensen): \(\log(\bar X) \neq \overline{\log X}\).

# La media del log NON è il log della media
log(mean(dati$WEIGHT))
#> [1] 4.977783
mean(log(dati$WEIGHT))
#> [1] 4.964609